Vortex waves in trapped Bose-Einstein condensates 
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We have theoretically studied vortex waves of Bose-Einstein condensates in elongated harmonic 
traps. Our focus is on the axisymmetric varicose waves and helical Kelvin waves of singly quantized 
vortex lines. Growth and decay dynamics of both types of vortex waves are discussed. We propose 
a method to experimentally create these vortex waves on demand. 
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I. INTRODUCTION 

Substance whirling around a line in space forms a vor- 
tex. Eddies in a flowing water and whirlwinds in the 
atmosphere are perhaps the most common examples of 
vortices encountered in nature, although such structures 
exist in nearly every imaginable medium including laser 
light and superconductors. Vortices are objects of fun- 
damental importance in fluid dynamics and they are in- 
herent to turbulent flows featuring normal mode excita- 
tions analogous to vibrating classical strings pQ. Lord 
Kelvin (Sir W. Thomson) derived a dispersion relation 
for a class of normal modes of classical hydrodynamic 
vortices [2]. A particular type of vortex excitations in 
which the vortex twists around its equilibrium position 
forming a helical structure have since become known as 
Kelvin waves [3 . Such chiral vortex vibrations are an ele- 
mental property of vortices ranging from tiny whirlpools 
in creeks to tornadoes of atmospheric scales. Another 
type of axially propagating excitation mode — coined a 
varicose wave — may exist in fluids confined to cylindri- 
cal geometry. These radially pulsating modes have the 
varicose shape, well known in the context of blood flow 
through veins in the cardiovascular systems [4]. 

Superfluids are distinguished from classical fluids by 
their ability of frictionless flow [5 j. Rotational mo- 
tion of superfluids is supported by quantized vortices 
[3] . Such quantum vortices possess many properties akin 
to their classical counterparts — including Kelvin waves 
IlllHiai^IinilllJIIlinilll. Superfluid Kelvin waves 
were first experimentally introduced in the context of 
liquid helium [3j [H [71 [9]. Later the discovery of Bose- 
Einstein condensates of dilute atomic gases [I5j [I6j [17] 
facilitated the creation and direct imaging of quantized 
vortices in highly controllable superfluid systems [I8j [19] . 
In particular, Kelvin waves on an isolated quantized vor- 
tex line were explicitly excited and imaged for the first 
time by Bretin et al. [12 . Their excitation method was 
described in terms of Beliaev decay of the quadrupole 
mode into a pair of Kelvin waves [T2J [2D] . Further in- 
sight for this process was recently obtained by direct dy- 
namical simulations [21]. Hydrodynamically it may be 
understood in terms of a parametric resonance between 
the quadrupole mode and a pair of Kelvin modes of oppo- 
site axial momenta, which induces an elliptical instability 
to the initially straight vortex line [2TJ [22] . The dynam- 



ics of superfluid Kelvin waves is crucial to understanding 
the details of quantum turbulence [23] [24] . It is also con- 
ceivable that the quantum gases would support varicose 
waves in which an axisymmetric pulsating perturbation 
propagates along the vortex line causing periodic modu- 
lation of the vortex core diameter. 

In this paper, we mainly focus on two types of vor- 
tex waves in trapped Bose-Einstein condensed gases. In 
the Sec. II we describe our numerical methods and in the 
Sec. Ill we first extend our previous results [21] on Kelvin 
waves followed by an introduction to varicose wave exci- 
tations. Finally, the Sec. IV is devoted to discussion. 

II. METHODS 

We consider a Bose-Einstein condensate composed of 
87 Rb particles trapped in a harmonic potential whose 
transverse and axial frequencies are uo x — uo y = co± = 
2tt x 98.5 Hz and uj z = 27r x 11.8 Hz, respectively, as in the 
experiment of Bretin et al. [12] . The number of particles 
in the condensate is 1.3 x 10 5 unless otherwise stated. The 
ground state and its dynamics are accurately described 
by the time-dependent Gross-Pitaevskii equation [25j [26] 

ihdt4>(r,t)= (H + V peTt (r,t) + g\^{r,t)\ 2 ) ^{r,t), (1) 

where the single-particle quantum harmonic oscillator 
operator is given by 

tt 11 V I / 2 2 , 22, 2 2\ /o\ 

= — 2^ 2 m ^ xX yV +C ^ Z )' ( 2 ) 

and the last term in Eq.Q accounts for the contact 
interaction of strength g between particles. The time- 
dependent potential V pert (r,t) may be engineered to res- 
onantly excite various normal modes by perturbing the 
ground state. Such collective excitation modes are also 
obtained within the linear response theory by solving 
the coupled Bogoliubov-de Gennes eigenvalue equations 
[25] [26] 

Cu q (r) ^ g^(r) 2 v q (r) = E q u q (r) 
£^(r)+^(r)*%(r) = -E q v q (r), (3) 

where u q and v q are the usual quasiparticle wavefunctions 
corresponding to the excitations with eigenenergies E q , 
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and the single-particle operator 

£ = tf + 2^(M)| 2 -/i (4) 

contains the condensate chemical potential \i. 

We solve the time-dependent Gross-Pitaevskii equa- 
tion by utilizing a parallel MPI implementation of a 
finite-element discrete variable representation method in 
a fully three dimensional cartesian coordinate system 
[27) [28] . For finding solutions to the Bogoliubov-de 
Gennes equations, we deploy the cylindrical symmetry 
of the axisymmetric single vortex stationary state. As a 
result we have a good azimuthal quantum number t thus 
reducing the effective computational effort by one spatial 
dimension. To diagonalize the matrix form of Eq.([3|, we 
use a hybrid basis consisting of radial Laguerre polyno- 
mials and finite-difference discretation in the axial direc- 
tion. This yields a matrix which has a computationally 
beneficial band shape. 



III. RESULTS 

By solving the eigenvalue problem for the axisymmet- 
ric single vortex initial state we obtain full information 
of the small-amplitude normal modes of our superfluid 
vortex state. Using dynamical simulations we can study 
the excitation and subsequent decay dynamics of vortex 
waves beyond linear stability analysis. Combining these 
two resources we are able to analyze excitation mecha- 
nisms, dispersion relations, and decay processes of the 
collective vortex excitations. In the following we sepa- 
rately focus on helical t = — 1 Kelvin waves and axisym- 
metric t = varicose waves. 



A. Helical £ — — 1 Kelvin waves 

Physical characterization. — Kelvin waves or their 
quanta — kelvons — are helical vortex waves propagating 
along the vortex line with an axial wave vector k 2 . The 
path so(C>£) traced by the core of a vortex executing an 
idealized classical Kelvin wave motion of an amplitude R 
may be curve parametrized as 

s o(C t) = R cos(k z ( + uot)e^ — R sm(k z ( + uot)e y + £e z , 

(5) 

where t denotes time and uo is the rotation frequency of 
the line element about the z-axis. A particle following 
such trajectory has a negative helicity. One may use the 
Biot-Savart law to obtain the velocity field [3] 

h f (s - r) x ds , a , 
V(r) = W |s -r|3 • ® 

The divergence in the integral may be removed by intro- 
ducing a cutoff ao = |r — so|, which causes the internal 
structure of the vortex core to be neglected. Within the 



local induction approximation the self-induced fluid ve- 
locity in the vicinity of the vortex core is given by [3] [29] , 

v loc (r)=Cs'xs", (7) 

where C is a system dependent constant and a prime 
denotes a derivative with respect to the curve parameter 
so that s" is the local curvature vector of the vortex and s' 
is tangent to the vortex line, pointing to the direction of 
the local vorticity vector. The Kelvin-Helmholz theorem 
guarantees the conservation of circulation which implies 
that a vortex tends to move with the fluid it is immersed 
in. Hence according to Eq.([7|) the local motion of the 
vortex line element is in the direction determined by the 
binormal of the vorticity and the curvature vectors. In 
a uniform system this then implies that the Kelvin wave 
vortex helix must rotate in the direction opposite to the 
fluid flow around the unperturbed vortex. 

However, in the presence of external fields such as a 
harmonic trapping field, the local curvature s" may be- 
come significantly modified. A particularly elusive exam- 
ple is a singly quantized vortex precessing about the trap 
axis in a harmonically trapped Bose-Einstein condensate 
[30] . In such system the vortex is bent even in the pan- 
cake limit and the effective curvature vector points away 
from the trap centre and therefore the local self-induced 
velocity, which determines the direction of vortex preces- 
sion, is in the same direction as the unperturbed con- 
densate flow. The Kelvin modes correspond to a small 
amplitude motion of the vortex core and the correspond- 
ing mode densities are highly localized within the core 
volume of the vortex. 

Dispersion relation. — The driving force of Kelvin wave 
motion involves velocity instead of acceleration and the 
resulting semiclassical long- wavelength, \a^k z \ <C 1, dis- 
persion relation uo{k z ) = hk 2 z In (l/|aofc 2 |) /2m depends 
logarithmically on the wave vector k z . In the case of 
an axisymmetric vortex line in a harmonically trapped 
Bose-Einstein condensate the kelvon dispersion relation 
is better described by [21] 

Hk 2 ( 1 \ 
uo(k + k z )=uo + ^-\nl yY\)' l rcfc *l <1 ' (8) 

where r c is a core parameter which is of the order of the 
vortex core size. The "continuum" branch described by 
the above formula is cut off at low frequency uoq < and 
long- wavelength fco ~ ^/D Z1 where D z is the Thomas- 
Fermi length of the condensate. In addition, the spec- 
trum contains low-energy bending modes which lie below 
the continuum branch. From a qualitative point of view 
modes with uo < rotate in the same direction as the 
unperturbed condensate flow, whereas modes with uo > 
rotate in the opposite direction. Modes with zero fre- 
quency remain stationary in the laboratory frame. Fig- 
ure [I] shows the computed kelvon dispersion relations 
for three different system sizes N = 1 x 10 4 (squares), 
iV = 5 x 10 4 (bullets), and iV = 5 x 10 5 (triangles). 
The solid curves are plotted using the dispersion rela- 
tion, Eq.([8|, where uoq and fco are determined from the 
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FIG. 1: Kelvon dispersion relation for different system sizes, 
N = lxlO 4 (squares), 5xl0 4 (bullets), and 5xl0 5 (triangles). 
The solid curves are plotted using the dispersion relation, Eq. 
J8|, with core parameter r c = 0.13 /im for all cases. The 
dashed line illustrates the sensitivity of the dispersion relation 
to the change of the core parameter as specified in the text. 

numerical results and the core parameter r c = 0.13 /mi 
for all cases. In contrast, the value of the healing length 
varies between 0.2 and 0.4 /mi. The dashed curve is plot- 
ted for N = 1 x 10 4 using r c = £ = 0.4 jam highlighting 
the failure of approximating the core parameter r c by the 
healing length £. However, all cases shown in Fig. [I] are 
still in the Thomas-Fermi regime and for very small sys- 
tems the dispersion relation shows much stronger finite- 
size effects. The modes below the continuum dispersion 
are the bending modes [21]. As the number of particles 
in the condensate is increased it becomes correspondingly 
longer and hence fco approaches zero. 

Excitation method. — An elegant way to excite Kelvin 
waves has been realized in the experiment of Bretin et 
at. [12 j. In this method, an elliptical rotating per- 
turbation first populates the quadrupole t = — 2 col- 
lective surface oscillation mode of the condensate thus 
reducing the orbital angular momentum of the system. 
Subsequently, the quadrupole mode population sponta- 
neously decays pairwise into t = — 1 Kelvin modes via the 
energy-momentum conserving Beliaev process in which a 
quadrupole quasiparticle interacts with the condensate, 
scattering to a pair of oppositely traveling Kelvin modes 
[12] [20] [21] • Hydrodynamically one may explain such 
process in terms of a parametric (subharmonic) reso- 
nance between a quadrupole mode and a pair of Kelvin 
modes, leading to an elliptical instability of the initially 
straight vortex line [2TJ [22] . 

The shortcoming of this resonant excitation method 
is that it is constrained to a specific kelvon. However, 
it is possible to circumvent this constraint by deploy- 
ing the core localized nature of the Kelvin modes [2T] . 
By applying an external pinning potential, such as pro- 
duced by a suitably detuned laser beam, at the vortex 
core thereby changing the local effective potential it is 




FIG. 2: Kelvin waves of different wave numbers excited using 
Vo/hv± = 6.0 (a) 4.0 (b), and -4.0 (c). The duration of the 
quadrupole perturbation is 35 ms for all cases. The frames 
(a)-(c) are for times 75, 86, and 112 ms, respectively. 

possible to control the energies of the kelvons with re- 
spect to the condensate ground state energy. Since this 
may be performed in a fashion which does not dramat- 
ically affect the resonant quadrupole mode frequency, it 
is possible to modify and choose the kelvon wave vector 
which resonates with the quadrupole mode. Quantita- 
tively such controlled kelvon excitation process may be 
achieved by using the time-dependent perturbation 

2 

Vp ert (r,*) = Vo^^eM-2r 2 /o-(z) 2 ) (9) 
+ emu; 2 L (cos(2nt)(x 2 - y 2 ) + 2sm(2Qt)xy)/2, 

where the first line models a Gaussian laser field and pro- 
vides the frequency shift to the kelvon dispersion. It is 
also worth noting that at finite temperatures the kelvon 
frequencies may experience an additional frequency shift 
[3Tj [32] . The weak potential, e = 0.025, in the second 
line of Eq.([9| may be realized by rotating an elliptically 
deformed trap and is responsible for the resonant exci- 
tation of the quadrupole mode. In Eq.([9|, Q = — 0.6cu_\_ 
is the rotati on frequenc y of the elliptical perturbation, 
cr(z) = do y/l + (z/zr) 2 and we have chosen cr = 2/mi 
with a Rayleigh range of = 20/im. It may be experi- 
mentally challenging to achieve focusing to such a narrow 
beam waist, however, the method also works for wider 
beams although the accuracy in determining the disper- 
sion relation may then suffer. Furthermore, it would 
probably be wiser to switch on the pinning beam already 
before creating the condensate in order to avoid super- 
fluous excitations. 

Frames (a)-(c) of figure [2] show a gallery of kelvons 
with different wave vectors which have been gener- 
ated by applying the above prescription with Vo = 
{6.0,4.0,-4.0} hcu_\_, respectively. The duration of 
the perturbation was set to 35 ms and the plots are 
for times 75, 86, and 112 ms after the beginning 
of the perturbation. The initial in-plane sinusoidal 
shape of the vortex s±(£, t) = 2R cos (k z () cos (cut) e x — 
2Rcos(k z () sm(cut)e y -\-(e z , seen in each frame of Fig. 
forms as a superposition of two helical waves traveling in 
opposite axial directions. It might also be possible to 
directly excite kelvons by using suitably tuned Laguerre- 
Gaussian beams or by applying an axially propagating 
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velocity field such as a frequency matched traveling op- 
tical lattice with an appropriately chosen lattice period 
which could provide a resonant coupling to Kelvin modes. 

Decay dynamics. — The decay of Kelvin waves may be 
studied by measuring the length of the vortex line 



S(t) = [ 

J-Z 



[dX(z,t)/dz] 2 + [dY(z,t)/dz] 2 dz 



(10) 

as a function of time t, where X(z,t) and Y(z,t) are 
the parametrized x and y positions of the vortex core, 
respectively. We define a reduced line length density 
L(t) = (S(t)/S(0)-i)mu±/h, where 5(0) is the reference 
length of a straight vortex. The quantity L{t)/L(t f ) is 
shown in Fig. (|3|, where S(t')/S(0) = 1.06 and the data 
is obtained for Vo/Huj± =0. At t « 0.1 s the vortex line 
shape rapidly changes from being straight to oscillating 
sinusoidally in plane. The increase in the line length may 
be estimated by an exponential function, Aex.-p(jujj_t) 
plotted in Fig. ^ using A = 0.001 and 7 = 0.3. Between 
t « 0.1 s and t ~ 0.3 s the line length diminishes as the 
kelvons decay to kelvons of smaller wave numbers. The 
final state is a bent vortex which causes the data in Fig. 
Q to saturate at a finite value. For our system, energy 
is conserved both during the growth and the decay of the 
line length. In the context of superfluid turbulence and 
Kelvin wave cascade, the decay of the vortex line length 
density is frequently modeled using a rate equation of the 
form [H [33 El EES M El] 



dL(t) 
dt 



(11) 



where, for a = 2, v has the units of an effective kinematic 
viscosity. This may be integrated to yield 



1 



Ht) = 

L(t') l + vL{t'){t-t')' 



(12) 



The decaying curve in Fig. Q is plotted using v/k = 0.2, 
where n = 27rh/m is the quantum of circulation. It is in- 
teresting to compare this result to, e.g., those in Refs. 
[34j [36j |37] . However, more detailed investigation is re- 
quired if one desires to draw a direct connection between 
the decay of Kelvin waves on a single vortex in a harmon- 
ically trapped Bose-Einstein condensate and the decay of 
vortex tangle in turbulent superfluid helium systems. 

Further insight to the excitation and decay of Kelvin 
waves may be obtained in terms of the multipole mo- 
ments of the system. Initially, the condensate wavefunc- 
tion attains strong quadrupole £ = — 2 moment since the 
elliptical rotating drive builds up population in the ap- 
propriate quadrupole mode. When the quadrupole mode 
population decays to £ = —1 kelvons, the wavefunction 
gains a corresponding dipole moment and finally, as the 
kelvons further decay to phonons and kelvons of longer 
wavelength, they reveal themselves in the monopole mo- 
ment of the wavefunction. Thus by analyzing the z- 
dependence of the real monopole £ = and dipole £ = 1 
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FIG. 3: Reduced length of the vortex as a function of time. 
Initially the vortex line remains straight until at 0.1 s its 
length rapidly increases due to the exponentially growing in- 
stability. The line length decay is modeled using the equation 
(12). 




FIG. 4: Fourier transformed dipole moment k\ (a) and 
monopole moment (b) of the system as a function of time. 
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y e mx,y,z)\ 2 dxdy (13) 



we may gain further information on the collective excita- 
tions present in the system. Figure |4ja) and (b) display 
the respective spatially Fourier transformed dipole £ = 1 
and monopole £ — signals 



kl(t)=f{Q e (z,t)-Q e (z,0)}, 



(14) 



where T denotes a spatial Fourier transform, obtained for 
Vo/fov± = 0. At t « 0.1 s a clear dipole moment signal 
emerges at k\ = 0.8 and 0.9 /im" 1 signifying the onset of 
the kelvon excitation, see Fig. |4ja). Soon afterwards, a 
strong monopole signal emerges at k° z = 0.2 and 0.4 /im 
due to the phonon emission from the decaying kelvons, 
see Fig. [4^b). Since the cylindrical symmetry is broken 
when the kelvons are excited, their presence in the system 
is manifest in both the monopole and dipole signals. 
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FIG. 5: Dispersion relation for varicose waves (boxes), sound 
waves (circles) and Kelvin waves (triangles). Solid (open) 
markers are plotted for even (odd) z-parity modes. The 
solid curves are the respective analytical varicon, phonon and 
kelvon dispersion relations discussed in the text. 



B. Axisymmetric £ — varicose waves 

Physical characterization. — We assign the term vari- 
cose wave — or varicon to denote their quanta — to de- 
scribe a class of axisymmetric, £ = 0, collective modes of 
the condensate. Here we focus on varicons of lowest fre- 
quencies which involve transverse motion and hence our 
varicose wave having zero axial momentum corresponds 
to the usual radial breathing (or monopole) mode [38 . 
The modes with finite axial wave number exhibit sim- 
ilar sinusoidal breathing oscillation with z— dependent 
phase. In these systems the collective modes are den- 
sity waves and therefore also the axial phonons produce 
similar density modulation along the vortex core due to 
the local density variations caused by the axially propa- 
gating sound waves [39] . In contrast to the kelvons which 
only exist in the presence of a vortex and are localized 
in the vortex core, varicons correspond to the transverse 
swell-squeeze motion of the condensate and exist also in 
the absence of the vortex. 

Dispersion relation. — Figure [5] shows the computed 
varicose mode frequencies (boxes) as a function of their 
wave vector. For comparison, we have also included both 
the phonon (circles) and kelvon (triangles) modes. The 
upper-most solid curve follows the frequency shifted free 
particle dispersion relation 



u(k z ) = uj + — — 



2m 



(15) 



describing a varicose mode in the presence of a vortex 
line. Here uoq corresponds to the fundamental monopole 
frequency of the condensate [381. The lowest curve is the 
kelvon dispersion relation, Eq.|8|, and the straight line 
is the linear phonon dispersion uo p = c\k z \, where c is the 
density averaged speed of sound. 
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FIG. 6: Varicose waves of different wave vectors k z = 0.2 (a) 
0.4 (b), and 0.6 /mi -1 (c). 
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FIG. 7: Wave vector as measured from the monopole moment 
according to Eq.(14) as a function of time. The varicon with 
k z = 0.4 /im _1 is first resonantly excited by the external 
perturbation. For later times, these varicons couple to the 
phonon branch, in particular to the one with k z = 1.1 /mi -1 . 
The inset shows the condensate at t = 48 ms. 



Excitation method. — Varicose waves may be reso- 
nantly excited using a perturbing potential of the form 

Vp ert (r, t) = emuo\{x 2 + y 2 ) cos(f)t) cos{k z z), (16) 

where Q estimates the varicon frequency and k z is the 
corresponding wave vector. In the limit k z — > the 
perturbation corresponds to a simple modulation of the 
strength of the radial trapping potential resulting in the 
excitation of the fundamental monopole breathing mode. 
Figure [6] shows a collection of varicose waves of differ- 
ent wave numbers excited using the potential given by 
Eq.([T6|, where the frequency and wavevector are chosen 
according to Fig. [5] The duration of the perturbation is 
24 ms and the frames (a)-(c) are plotted for times 48, 26, 
and 27 ms, respectively. 

Decay dynamics. — Figure [7] displays the wave vectors 
of the varicose waves as measured from the monopole 
moment Q° according to Eq.([l4|, as a function of time. 
The resonant external perturbation first excites the vari- 
con with k z = 0.4 /im -1 which then couples to phonons 
and in particular to the one with k z = 1.1 /im -1 , see 
also Fig|7l The varicose wave perturbations propagate 
in both directions along the vortex axis. On arriving 
at the low-density ends of the condensate, the varicose 
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waves produce ring currents which become pinched off 
from the mother condensate as seen in the inset of Figj7[ 

IV. DISCUSSION 

In conclusion, we have studied axially propagating col- 
lective excitations of quantized vortex lines in harmon- 
ically trapped elongated Bose-Einstein condensates. In 
particular, we have focused on two types of axial vortex 
waves: helical Kelvin waves and axisymmetric varicose 
waves. We have studied the excitation mechanisms, dis- 
persion relations and decay processes of both types of 
collective excitation modes. We have employed efficient 
computational methods to solve: (i) the time-dependent 
Gross-Pitaevskii equation for the full three-dimensional 
quantum dynamics of the vortex motion and (ii) the 
Bogoliubov-de Gennes equations for the stationary ax- 
isymmetric single vortex eigenstates. The combination 
of these two tools provides us with a comprehensive de- 
scription for this zero temperature system. 

The Kelvin wave dispersion relation at the long wave 
length limit has both positive and negative energy eigen- 
modes with respect to the condensate energy and these 
have distinct physical manifestations: positive modes 
counter-rotate with respect to the unperturbed superflow 
while the negative modes are co-rotating. The kelvon fre- 
quencies may be engineered using external pinning poten- 
tials. In particular, it is possible to choose the sense of 
rotation of particular kelvon and to excite kelvons of spe- 
cific wave numbers. Experimentally, Kelvin waves may 
be created indirectly by exciting the quadrupole mode 
of the system and allowing it to resonate with a pair of 



kelvons. 

The only dissipative mechanism operational in the zero 
temperature limit which may return a turbulent super- 
fluid system to a laminar flow state is a reduction of 
the total length of the vortex lines. It is believed that 
in an isolated vortex the Kelvin wave cascade process 
is responsible for this transfer of the excess energy from 
small to large wave numbers and eventually phonon ra- 
diation is expected to absorb the surplus energy of the 
turbulence [231 EH EQl HU H21 [43]. In contrast to this 
scenario — if applicable to our freely decaying single vor- 
tex case — here the cascade proceeds predominantly to- 
ward smaller wave numbers, i.e., toward longer wave- 
lengths and smaller kelvon frequencies. 

In addition, we have studied axisymmetric excitations 
which exhibit a varicose shape. The long wavelength 
dispersion relation of the varicose waves studied here is 
quadratic in the wave vector. The varicose waves are 
readily visualized as periodic modulations of the vortex 
core diameter along the vortex axis, although the modes 
themselves exist independent of the presence of a vortex. 
The varicons are found to decay via coupling to phonon 
modes of higher wave vector. Detailed understanding 
of quantized vortex waves and their dynamics may shed 
light on open questions regarding various vortical systems 
including turbulent superfluids. 
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